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CALIFORNIA 


ABSTRACT 


The finite element method is used to study the structural 
response of a submarine to an acoustic blast wave. Encounter 
geometry is restricted to the case of a plane wave front 
parallel to the axis of a right-circular cylindric structure 
It is shown that the blast wave propagation may be studied 
separately from the structural response. For the case 
considered this separation allows a two-dimensional formula¬ 
tion of the propagation problem v/hile retaining a three- 
dimensional structural model. Propagation results are found 
to agree well with analytic results obtained by others. 
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I. INTRODUCTION 


Consider a nonrigid structure submerged in a fluid. 
Pressure in the fluid may induce response motions of the 
structure. Conversely, motion of the structure affects 
the pressure field. Such a system is called a coupled 
elasto-hydrodynarnic system. 

Many practical engineering problems belong to this 
/ 

category. A resei’volr-dam system subjected to earthquake 
motion, oscillations in rocket fuel systems, a vibrating 
ship’s hull, and structures subjected to blast loadings due 
to explosions are familiar elasto-hydrodynarnic problems. 

The particular elasto-hydrodynarnic system under consid¬ 
eration in this paper is a submarine subjected to a blast 
wave resulting from an underwater explosion. The submarine 
is modeled as a ring-reinforced, neutrally buoyant, circular 
cylinder with rigid bulkheads. It is assumed that the blast 
wave may be modeled as a plane acoustic wave with wave-front 
parallel to the longitudinal axis of the submarine. 

The interaction of an acoustic v;ave with an obstacle has 
been the subject of many investigations. Sette [1] first 
computed the pressure distribution on a rigid cylindrical 
surface due to a step-pulse. The response of a cylindrical 
shell to a plane step (or decaying) shock wave has been 
treated by various authors [2-5]. 

Due to the complex nature of the coupled response a 
realistic analytic formulation of this problem is intractable. 
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However, numerical methods, particularly the finite element 
method, provide a satisfactory means of solution. 

The finite element formulation of structure-fluid 
dynamics was first introduced by Zienkiewicz, Irons an^ 

Nath [6]. They demonstrated that the coupled natural 
frequencies and mode shapes of systems could be accurately 
determined by this method. 

It is the purpose of this paper to demonstrate that the 
finite element method, as forwarded by Zienkiewicz and 
Newton [7], is applicable to the solution of the interaction 
of a blast wave with a rigid structure. Using the finite 
element formulation, a superposition theorem is developed 
which affords a method of applying the rigid cylinder results 
to a realistic, but simple, three-dimensional model of the 
submarine. 
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II. THE FINITE ELEMENT METHOD OF SOLUTION 


A. THE STRUCTURE DISCRETIZATION 

The finite element discretization of the structure 
problem is well known and is described in Zienklewicz's 

S 

text [8]. The discretization process results in an 
assembled system of equations of the form 

V 

CM]{6} + CC]{6} + [K]{6} = {R} . (1) 

In the above equation [M], [C] and [K] are respectively the 
mass, damping, and stiffness matrices of the structure 
calculated in the proper, consistent, manner, {6} is a vector 
of nodal displacements, the dots indicate time differentiation, 
and {R} is a vector of generalized nodal loads. 

It is convenient to divide the generalized nodal forces 
into two parts (after Zienklewicz and Newton [?]) such that 

{R} = {F} + {P} . (2) 

The vector {F} is due to external forces and {P} is due to 
the fluid pressure on the interface. For the i^^ structural 
node 


P = / N ' p dS , (3) 

1 g 1 

where N^^' is the appropriate shape function defining the 
displacement pattern in the direction normal to the boundary, 
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p is the pressure on the structure-fluid interface S, and 
the integration is over the Interface. 

B. THE FLUID DISCRETIZATION 

From first principles it can be shown [6,9] that the 
wave equation 

V^P = 4r p , (H) 

c 

2 

where V is the Laplace operator, p is the hydrodynamic 
pressure in excess of static pressure and c the acoustic 
velocity of propagation, together with the necessary 
boundary and initial conditions are the governing relations 
for the hydrodynamic response of the fluid. 

The spatial discretization of the fluid problem is 
accomplished by a direct application of the Galerkin weighted 
residual process . If at any instant of time p is approximated 
by 


p - N.p. (summation on J = l,...,m) , (5) 

J J 

where m is the number of nodes, N. the j shape function 

J 

chosen to describe the spatial variation of p over the region 

and p. is a set of nodal pressure values which are time 
J 

dependent, then the 1^^ v/elghted residual equation is 

/ N. [V^(N.p,) - ^ (N.p,)] dR = 0 , (6) 

R 1 JO J J 
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where is the shape (weighting) function [8, Ch. 3] of 
the Galerkin process and R is the fluid region under 
consideration. Employing Green's theorem, the v/eighted 
residual equations may be transformed to 


/ 

R 


8x 


8N. 

•r—^ + 
8x 


8N. 8N 

_1 

3y 8y 


3N. 8N. 
_ 1 ,1 

dz dZ 


dR p. 


2 


/ N.N 

R 1 J 


dR p. 


aN. 

/ N. —^ 

S ^ 


dS p.. 
3n 


= 0 


(7) 


This system of equations can be written in matrix form as 

[Q]{p} + [H]{p} = {B} . (8) 

The matrices [Q] and [K] are the assembled "inertia" and 
"stiffness" matrices of the fluid. They are defined on the 
element level by: 




/ N N dR 

T? J ^ 


(9) 



/ 


R, 


9N. aN. aN. aN, aN, aN. 

_^ ^ ^ 1 ^ ^ i 1 

ax ax ay ay az az e 


( 10 ) 


where R^ denotes the element region. 

The vector {B} results from the surface integral in (7) 
and can be v;ritten for the 1^^ surface node as 


b 


1 


/ 

S 


M. 


an 


dS 


( 11 ) 
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where n Is a local coordinate in the direction of the 
outward normal. By means of this vector it is possible to 
apply the boundary condition on the normal pressure gradient, 
i .e., 


9P 

9n 




( 12 ) 


where is the fluid particle velocity normal to the 
boundary and p is the fluid density. It is noted that for 
a fixed boundary “ 0 the corresponding b^ s 0 . 

A ’nonreflecting* or radiation boundary condition, 
necessary to model numerically an infinite fluid region, vjas 
developed by Zienklev^lcz and Newton [7]- They demonstrated 
that for a plane wave normally incident upon a plane boundary 
the condition to be satisfied by p on the boundary S is 



The corresponding boundary Integral (11) is 


(13) 


I ^ % H os = 


i / N^Nj dS Pj 


( 1 ^) 


This leads to an additional "damping” term [D]{p} in the 
matrix equation (8), i.e.. 


[Q]{p} + [D]{p} + [H]{p} = {B} 


(15) 
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[D] is the assembled fluid "damping" matrix defined on the 
element level as 



( 16 ) 


where denotes the external boundary of the element. 

C. THE COUPLED STRUCTURE-FLUID PROBLEM 

The vectors {P} of the force term of equation (2) and 
{B} of equation (15) determine the coupling of the fluid 
and structure in the discretized problem. Recall from 
equation (3) that the nodal forces on the structure due to 
the pressure are 

P. = / N.' p dS = / N.'N.p. dS = [ / N.'N. dSj p. 


s 


s 


or 


(17) 


{P} = [L] {p} 


where [L] is the assembled matrix of element contributions 



(18) 


From equation (12) 




















The forcing term (b) In equation (15) becomes 





3 


or 


{B} = - [S]{6} 


(20) 


with [S] the assembled matrix of element contributions 


" p / NlNj dS^ • (21) 

Observe that 

p[L]'^ = [S] . (22) 

The complete elasto-hydrodynamic system has been 
formulated and can be summarized. Rewriting equation (1) 
with (2) and (17), 

CM]{6} + [C]{6} + [K]{6} = [L]{p} + {F} (23) 

describes the structure behavior. Similarly, equation (15) 
v/ith (20) and (22) can be written as 

CQ]{p} + [D]{p} + [H]{p} = - p [L]^{6} . 


This governs the fluid response. 
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The complete system of equations can be written as 



Using general finite element techniques It Is possible 
to generate and subsequently solve this unsymmetrlcal system 
of equations for the submarine-fluid problem under considera¬ 
tion. However, an alternative method based on this general 
formulation V7lll be utilized to effect a solution using the 
following superposition theorem. 

D. SUPERPOSITION THEOREM 

It will be shov;n that the coupled elasto-hydrodynamlc 
problem formulated above may be solved In two stages and that 
the structure response thus found Is the same as In the 
complete unseparated problem. In the first stage the structure 
Is constrained against responding to the loading Imposed by 
the Incident pressure pulse. In the second stage the structure 
and the fluid are Initially at rest. External loads are 
applied to the structure to simulate the loading by the 
pressure pulse of the first stage and the structure and fluid 
responses are found. 

1. The Complete Problem 

Equations (23) and (2^) govern the structure and 
fluid response subject to the Initial conditions 
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{<5(0)} 

= {^(0)} = 0 , 

(26) 

{p(0)} 

= {a} , {p(0)} = {b} , 

(27) 


where {a} and {b} are given vectors. It is assumed that 

{P} = 0 . ( 28 ) 

The solution of the complete problem requires that equations 
( 23 ) and (2^) with equation (28) be solved for {6(t)} and 
{p(t)}, subject to the Initial conditions, equations (26) 
and (27). 

2. Stage On e 

The dependent variables will be {p'“\t)} and i6^^^(t)}. 
Recall that In the first stage the structure Is constrained 
against motion, l.e., 

{6^^^t)} = 0 . (29) 

For the fluid, equation (2^) becomes 

[Q]{p^^h + [D]{p^^h + [H]{p^^h = 0 , ( 30 ) 

with the Initial conditions 

{p^^^(O)} = {a} , (p^^^O)} = (b) . (31) 
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Equation (23)> governing the structure response, will be 
satisfied by (29) If It Is required that 

{F^D} = _ [L]{p^^h . (32) 

For the first stage one need only solve equation (30) 
for {p^^^(t)} subject to the Initial conditions of equation 
(31). Is subsequently found from equation ( 32 ). 

3. Stage Two 

The system Is Initially at rest and the structure Is 
externally loaded to simulate loading by the pressure pulse. 
The dependent variables for this stage are {p “^(t)} and 
{6^^^(t)}. For the fluid 

[Q]{p^^h + [D]{p^^h + [H]{p^^h = , (33) 

with the Initial conditions 

{p^^^O)} = {p^^^O)} = 0 . (3^) 

For the structure 

+ [C]{(S^^h + = [L]{p^^h + {F^^h, (35) 

with the Initial conditions 

{6^^^(0)} = {6^^^(0)} = 0 . (36) 
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specify further that 


{F^2)j ^ ^ (27) 

For the second stage of the solution it is necessary 
to solve equations (33) and (35) for {p*‘^^(t)} and {6^^^(t)} 
with equation (37) subject to the Initial conditions of 
equations (3^) and (36). 

. Superposition 

It is now asserted that 

{p(t)) = {p'^>(t)> + {p<^'(t)) , (38) 

and 

{6(t)} = {6^^^(t)} + {6^2^(t)} . (39) 

In virtue of the conditions imposed on the first 
stage pressures {p^^^(t)} and displacements {6^^^(t)} , 
equations (29) - (32), and on the second stage pressures 
tp^^^(t)} and displacements {6^^^(t)} , equations (33) - ( 37 ), 
it follows that {p(t)} and {6(t)} as given by equations (38) 
and (39) satisfy the governing equations for the complete 
problem, equations (25) - (28). This establishes the super¬ 
position theorem for the coupled elasto-hydrodynamic problem. 

Implementation of the superposition theorem does not 
require that the rigid body pressure distribution or the 
structure response be obtained by the finite element method. 
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It Is noted that for the submarine structural model under 


consideration a three-dimensional analysis is necessary. 
However, the superposition theorem permits the rigid body 
pressure distribution to be obtained from a two-dimensional 
analysis. 
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III. WAVE PROPAGATION BY THE FINITE ELEMENT METHOD 


To evaluate numerically the Interaction of a wave with 
a structure it is necessary to propagate the blast wave 
through the contiguous fluid region. Before proceeding with 
the complicated two-dimensional interaction of a blast wave 
with a rigid cylindrical structure, the characteristics of 
the one-dimensional wave propagation are studied. The 
techniques developed and utilized for the one-dimensional 
case are subsequently employed in the two-dimensional study. 

A. ONE-DIMENSIONAL V/AVE PROPAGATION 

The numerical results obtained by the finite element 
method for the one-dimensional vjave propagation problem are 
easily compared to well known theoretical values and thus 
provide insight into the behavioral characteristics of the 
techniques employed. 

The physical situation to be considered is a semi-infinite 
fluid strip of unit cross-sectional area. Through this region 
a pressure wave of known shape is propagating with acoustic 
velocity c. A rigid structure is located at x = 0 and the 
radiation boundary condition is applied at x = L. Equation 
(2h), with {6(t)} = 0, 

[Q]{p} + [D]{p} + [H]{p} - 0 , (^0) 

governs the fluid response for the case. The Initial 
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conditions {p(0)} and {p(0)} are determined by the wave form 
under consideration. 

1. Blast Wave Model 

A blast wave, which is generated by an underwater 
explosion, is characterized by an extremely rapid rise in 
pressure to a maximum value p*, followed by decreasing 
pressure. For the purpose of investigating blast wave 
propagation it would appear that a step pulse of finite 
duration v;ould be a satisfactory wave model. However, it is 
evident that to model a moving discontinuous function by a 
finite fixed spatial discretization is not possible. Conse¬ 
quently, for the purpose of this investigation the blast wave 
is modeled as a ramp rise to a constant value, i.e., a ramp- 
step wave. Although the ramp-step wave is of primary interest, 
this investigation is not limited exclusively to that wave 
form and illustrative results of tests made using different 
wave forms are also given later. 

2. Element Selection 

To study a significant portion of a blast wave, which 
is propagating at the acoustic velocity, requires a spatial 
domain of considerable extent. This is necessary to avoid 
spurious effects caused by reflections from the (artificial) 
boundaries. Accordingly, the required linear dimensions of 
the region are approximately proportional to the time interval 
during v;hich the structure-fluid response is sought. It is 
evident that, if a fine subdivision of the one-dimensional 
fluid strip is required to represent the v/ave form under 


25 





consideration, the number of degrees of freedom in equation 
(i<0) may become prohibitively large when the same subdivision 
is applied to the two-dimensional region. 

Two types of finite elements, linear and cubic, are 
considered for the spatial discretization of equation (^0). 

As their names imply they respectively provide for a linear 
and cubic variation of the pressure within an element. When 
propagating smooth wave forms represented by a significant 
number of elements, e.g., p(x) = p* sin (ux/L) or 
p(x) = p*cos (irx/L), both the linear and cubic elements 
produce excellent numerical results. However, In choosing 
the type of finite element best suited to represent spatially 
a propagating blast wave it is necessary to consider: a) the 
general shape of the wave; and b) two instantaneous finite 
element representations of the wave front as it propagates 
across a region. 

As established above, the general shape of the vmve 
of primary Interest is the ramp-step wave. For the second 
consideration, assume that the node spacing s^ is equal for 
both element types. Due to the character of the linear and 
cubic elements note that three linear elements are required 
for each cubic element in the spatial representation of the 
wave front. Nov; consider a ramp wave front that at one 
Instant in time is represented spatially by four nodal values, 
l.e., the ramp rise spans one cubic element. At this 
particular Instant both element types are capable of an exact 
representation of the ramp v;ave front. For the second Instant 
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in time, let the wave advance a distance equal to one node 
space. At this instant of time the linear elements continue 
to provide an exact representation of the wave front. 

However, when the wave is in this spatial position relative 
to the cubic elements, they are required to fit a cubic 
polynomial through a function which has a discontinuous first 
derivative. This results in a distortion of the v;ave front. 

The above qualitative discussion explains the 
observed superior performance of the linear element when 
used in the spatial discretization of equation (^0). It also 
serves as a reminder that higher order finite elements per 
se do not guarantee superior results and may indeed lead to 
poorer results. 

3• Time Integration 

The time Integration of the governing equation (40) 
is computed using two finite differencing techniques, the 
Nev/mark-3 method as adapted by Chan, Cox and Benfield [10] 
and the Houbolt method [11]. Both provide a numerical method 
of determ.lning the time dependent response of a differential 
equation of the form 

mw + dw + kv; = f (t) , (^ 1) 

where m, d and k are constants. The difference equation 
proposed by Chan, et al. [10] is 

(m + I d + eh^k)v;j^ = (2m’- (l- 26 )h^k)v/^_j^ 

- (m - I d + eh 2 k)w ^_2 ^ Vl ^ V 2 ^ ’ 
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where h is the time step and 6 is a parameter which may have 
any value from 0 to 1/H. The value B = 1/12 is chosen for 
calculations using the stability requirements established by 
Newmark [12]. 

The starting equation is 


3 

(m + I d + 3h^)w^ = (m + I d - (| - B)h^k - - B) Wq 

2 2 

+ (m - hWg + 

+ ( (^- B) + (I- - B) ) h\ . (43) 


where and Wq are initial values. 

The Houbolt method [11] is developed from a consider¬ 
ation of a cubic curve that passes through four successive 
ordinates to obtain the difference equation 

(2m + ^ dh + kh^)w^ = (5m + 3dh)w^^^ - (i|m + | dh)w ^_2 

+ (m + ^ dh)w ^_2 + h^f^ . 

The starting procedure for the Houbolt method employed 
in the one-dimensional propagation problem is as follovrs: 
assume that the wave form is known at t = 0 and that the wave 
is propagating undeformed through the fluid region vjlth the 
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acoustic velocity c. The pressure distribution p can be 
described as 

p(x,t) = f(x + ct) (I 15 ) 

for a wave moving in the negative x-direction. It is apparent 
that once the wave form is specified at t = 0 as a function 
of X the required preceding values of pressure may be calcu¬ 
lated directly from equation (45) and thus provide exact 
starting values for equation (44). 

It is noted that the Houbolt method is theoretically 
stable for all step sizes, while the Newmark-3 method becomes 
unstable if the limiting step size established by Newmark [12] 
is exceeded. 

As an example of the results obtained by the methods 
outlined above consider a fluid strip in which the acoustic 
velocity c = 5000 ft/sec. The region of interest extends 
5000 ft. from a rigid structure located at x = 0 and is 
divided by 31 equally spaced nodes. The integration time 
step h = 0.01 sec. At time t = 0 the wave form has the shape 
and position in the region shown in Pig. 1. As the wave 
strikes the rigid structure at x = 0 it is reflected and 
propagates to the right where the radiation boundary condition 
is applied. 

After 100 time steps (t = 1.0 sec.) the wave V 7 ill be 
positioned as shown in Pig. 2. The corresponding numei’ical 
results are shown for the Newmark-3 integration with cubic 
elements and the Houbolt integration - linear element 
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Fig. 1. Pressure-vs.-distance at t = 0 
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combination. The other combinations of integration techniques 
and elements give essentially identical results. 

It is evident that to determine the optimum integra¬ 
tion technique it is necessary to establish criteria to judge 
the results of tests. The criteria chosen are: a) the 
pressure-time history at the structure-fluid interface 
(impulse) and b) retention of wave shape after propagation 
through the fluid. 

If the rise time of the ramp is chosen as the 
parameter necessary to describe the wave front, criteria for 
selecting the time h and the node spacing s^^ can be developed. 

Assume that to describe the slope of the ramp rise 
at least tvjo linear elements are required, i.e., 

2 s 

T > —. (i<6) 

r — c 

Intuitively this restriction on the node spacing seems 
reasonable because it will provide that at least two spatial 
locations in the fluid strip are used to describe the slope 
of the ramp. 

It also seems reasonable to restrict the time step 
in such a way that the wave front cannot cross an element in 
one•step, i.e ., 

h < — . (^<7) 
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The above restrictions are proposed as a general 


guide in the initial selection of the nodal spacing and the 
integration time step. They are not intended to establish 
restrictions which cannot be violated. As shown belovj, 
satisfactory results may be obtained when the restriction on 
(^6) is violated. 


Using the criteria established above to judge the 


effectiveness of the integration technique, numerous tests 
were made to determine which time Integration technique is 
better suited for the purpose of propagating a ramp-step wave 
Consider, for example, a triangular shape wave that 
has the following characteristics: 



s 


c 


T 


r 



where is the time required for the wave to decay from the 
maximum value p* to zero. The pressure-time histories at the 
rigid structure-fluid interface obtained by both integration 
methods are shown in Pig. 3- The result shown in Fig. ^ is 
obtained by employing a cubic finite element discretization 
with the Houbolt integration technique for the same wave form 
Consider a ramp-step wave with the characteristics: 
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Fig. 3. Pressure-vs.-Time at the Rigid Structure-Fluid Interface 
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Fig. 4. Pressure-vs.-Time at the Rigid Structure-Fluid Interface 

































Fig. 5. Pressure-vs.-Tine at the Rigid Structure- 
Fluid Interface (rarnp-step v:ave). 
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The resulting pressure-time histories at the fluid-structure 
Interface are shown in Fig. 5. From this and other results 
it is concluded that the Houbolt integration technique is 
better suited for the propagation of the blast wave through 
a fluid region. 

B. TV/O-DIMENSIONAL WAVE PROPAGATION 

Based upon the results of the one-dimensional study, the 
linear finite element and the Houbolt integration technique 
are selected to evaluate numerically the pressure-time 
history at the rigid cylindrical structure-fluid interface 
resulting from a ramp-step wave. The geometry of the problem 
is shown in Fig. 6. 




0 


L 


Fig. 6. Two-dimensional geometry 
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The rigid cylindrical structure, of radius a, is centered 
at the origin of the x,y coordinate system. The diagram- 
matically shown ramp-step wave is propagating through the 
fluid in the negative x-direction with the acoustic velocity 
c. Due to symmetry only the upper half plane is shown and 
the rigid body boundary condition is applied along the 
x-axis. The additional necessary boundary and initial 
conditions are specified below. It is convenient to define 
a nondimensional time parameter t = ct/a, where t = 0 is the 
time of incidence of the wave at r = a, 0 = 0°. 

1. Boundary Conditions 

The radiation boundary condition as developed earlier 
is valid only for a plane wave normally Incident upon a 
plane boundary. Clearly the v^ave reflected from the rigid 
cylinder will not be normally incident upon the far fluid 
boundaries. It is necessary that these boundaries be 
located a sufficiently large distance from the structure to 
Insure that the results at the Interface are not affected. 
This investigation and Sette's [1] results both indicate that 
the transient response of the blast wave-cylindric structure 
Interaction is essentially completed in 8a/c time units after 
incidences at r = a, 0 = 0°. As a consequence the far fluid 
boundaries are located sufficiently far from the structure 
to prevent any boundary reflected \javes from reaching the 
structure in the time specified. For programming convenience 
the rigid body boundary condition is applied to the boundary 
at y = L and the radiation boundary condition is applied at 

X = ±L. 
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A significant characteristic of the radiation boundary 
condition v;as observed during the two-dimensional study. It 
v/as observed that if the constant value portion of the ramp- 
step wave is applied to the boundary nodes at x = L for three 
successive time increments, then the radiation boundary 
condition continued to let the wave enter the fluid region of 
interest. In addition the interaction of the v;ave reflected \ 
from the structure with the boundary does not affect the 
incoming wave. This observation permitted a significant ’ 

reduction of the fluid region under consideration, resulting 
in a proportionate reduction in the number of degrees-of- 
freedom required for the spatial discretization of equation (40). 

2. Two-Dimensional Mesh Consideration 

In considering the two-dimensional spatial representa¬ 
tion (mesh) of the fluid field by the finite element method, 
it is necessary to Incorporate both the physical geometry and 
the vrave propagation aspects of the problem. As described in 
the one-dimensional study the integration time step is 
logically limited by the node spacing. Equation (4?) governs 
this dependence. The physical dimension of the cylindric 
structure, v/hen considered in the context of the advancing 
wave front, also restricts the nodal spacing. Clearly, if 
the node spacing s^ and structure radius "a" are chosen to 
be of the same order of magnitude and the upper limit on the 
time step is utilized, the wave front could advance across 
the structure in a few time steps and the pressure-time 
history at the interface v;ould be of suspect value. It is 
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therefore necessary to produce a mesh which affords a degree 
of fineness near the structure compatible with the integration 
time step and wave characteristics. Additionally, computer 
storage and processing time restrictions must be considered. 

After extensive experimentation the mesh shown in 
Pig. 7 was developed. Due to geometric symmetry about the 
y-axls only the region of the first quadrant is shown. This 
mesh has the following notable characteristics: 

a) Elements which form the cylindrical structure boundary 
are ’special' 5 noded, linear-parabolic elements. At 
the interface the 'special' elements represent the 
rigid boundary by a parabolic arc and provide for a 
transition to the linear elements in the fluid field. 

The interface nodes are equally spaced at angular 
Intervals of 7.5°. 

b) The 'near fluid' region is characterized by a radial 
node spacing equal to the 'far fluid' node spacing and 
an angular node spacing of 15°• 

c) The transition from the cylindrical 'near fluid' region 
to the rectangular 'far fluid' region is accomplished 
with a minimum number of distorted and degenerate 
elements. 

d) The uniform rectangular 'far fluid' region is charac¬ 
terized by node spacing s^ = — a . The appropriate 
boundary conditions as discussed above are applied at 
X = ±6a and at y = 6a. 
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The overall characteristic of the mesh is element 
size uniformity. Results indicate that such uniformity in 
both spatial directions is desirable for satisfactory two- 
dimensional wave propagation. 

3. Integration Starting Procedure 

The Houbolt integration starting procedure, equation 
(^15), used for the one-dimensional case assumed a general 
solution of the vjave equation (^) of the form 

p(x,t) = f (x + ct) + g(x-ct) (^8) 


with g(x-ct) = 0. This solution is valid assuming that 
f'(x) and f"(x) exist. However, the ramp-step wave has 
a discontinuous first derivative and singular higher 
derivatives. To avoid this inconsistency assume that the 
ramp wave front takes the form of a cubic curve at t = 0. 


P 


^ + ax' + 3(X')^ 


(ii9) 


with = a + 33(x')^ 

d X 


(50) 


where x' = x - Xq . 

If p = p* and 1^ = 0 at x' = C as shown in Fig. 8 

d X 

it is then possible to solve equations (^19) and (50) for the 


unknovm coefficients 

























































Pig. 8. Cubic Wave Front 


a = 


3p* 

TT 


and 


® ■ 1,(3 • 


or 


P = o P’ 


n + 3 2L!. 


- (—)^ ) 


(52) 


By substituting into (4 8) and differentiating v;ith respect 
to time 
P 


t=0 


= 3 1^0 ( 1 - (^)=) . 
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(53) 







































Equations (52) and (53) determine the initial conditions 
necessary for a proper formulation of the boundary value 
problem in the region of the cubic ramp wave front. 

To Incorporate the initial conditions into a starting 
procedure for the Houbolt integration technique consider a 
Taylor series expansion about t = 0 

P “ Pq + Pot + I Pot^ . (S'!) 


where p^ is p 


t=0 

If h is the time step 


h2 •• 

p_l = pp - hpp + -^ Pq , 


• h^ ” 

P+1 = Pq + hpQ + — Pq . 




The vector (Pq) is deduced from equation (^0), l.e.. 


(Pq) = (-[C]{Po) - [H]{Po)) • (56) 


Thus, equations (55) with equations (56), (52) and (53) 
provide a consistent starting procedure for the two- 
dimensional propagation problem employing the Houbolt 
integration technique. 

To minimize errors in the numerical solution of the 
pressure-time history, the \mve front is positioned in the 
fluid so that at t = 0 the cubic-ramp has its initial (leftmost) 
point at r = a, 9=0°. 




. Results 


Consider the following integration, wave and mesh 
characteristics: 


h = 


_ a 




T = 0»97a 
^r c ^ 


S = 2 ^ 


Implementation of boundary conditions, mesh and 
starting procedures described above results in the pressure- 
time histories at 6 = 0° and 6 = l80°, shown in Pigs. 9 and 
10. The solid line represents Sette’s [1] results displaced 
along the time axis to coincide with the incidence of the 
constant portion of the cubic ramp-step wave. Additional 
curves are given for 6 = 30°, 60 °, 90°, 120°, 150° in 
Appendix A. No comparison is available for the curves given 
for 120° and 150°. Tabulated values of pressure vs. time at 
each structural node are given in Appendix B. Figures 11 and 

12 provide a comparison of results when the time step is 

a 0 80 a 

increased to h = and the rise time decreased to t = —- 

3c r c 

The results shown in Figs. 9 and 10 and Appendix A are 
in good agreement with those calculated by Sette [1]. This 
agreement establishes that the Houbolt [11] integration 
technique and the linear finite element are a suitable 
combination for evaluating numerically the interaction of a 
blast vjave v.'lth a rigid cylinder. 

















































Fls- ?. Fressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (6 = 0°) 
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Fig. 10. Pressure-vs.“Time at the Rigid Cylinder- 
Fluid Interface (0 = l80°) 
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Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (0 - 0°) 


Fig. 11. 
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Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (0 = l80°) 
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Pig. 12. 

























IV. THE SUBMARINE STRUCTURAL MODEL 


A simple three degree-of-freedom submarine structural 
model Is developed. This model demonstrates the implementa¬ 
tion of the superposition theorem using the pressure-time 
history obtained from the two-dimensional finite element 
solution of the Interaction of a blast v;ave with a rigid 
cylinder. 

A single, ring reinforced section of a subm.arine, 
terminated at both ends by bulkheads, is modeled. The effects 
of the stiffening rings are included by treating the shell 
as orthotropic. Each bulkhead is rigid in its ov;n plane. 

The added mass effect due to the surrounding water is 
included by considering the water to be incompressible and 
using strip theory. 

The three modes considered are: 1) rigid body mode, 

2 ) shear beam mode and 3) buckling mode. The equations of 
motion take the form 

[M]{q} + [K]{q} = {f} , (57) 

where [M] and [K] are the structural mass and stiffness 
matrices respectively, {q} the vector of generalized displace¬ 
ments and {f} the vector of generalized forces. 

A point on the shell. Fig. 13, is located by specifying 
z and 6. 
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Pig. 13 . Submarine Structure Geometry 


The following symbols are used to describe the structure 
model: 

a = shell radius, 

h = shell thickness 

s 

= bulkhead spacing, 
u = radial displacement 

V = tangential displacement 

w = axial displacement 

E = Young's modulus of elasticity, 

G = shear modulus of elasticity 

V = Poisson’s ratio. 

Values of the structural parameters used in calculations 
are specified in Appendix C. 
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The generalized modal displacements are defined as: 

1) Rigid body mode, 

u = -q^ cos 6 , V = sin 6 , w = 0 . (58) 


2 ) Shear beam mode, 

TT Z 

u = -q 2 sin — cos 6 , 
w = 0 . 

3) Buckling mode, 

. 2 TTZ 

u = -q^ sin -j- cos 26 
w = X q^ sin cos 26 


TTZ . - 

V = q 2 sin -j- sin 6 


(59) 


V = 


1 . 2 
2 sin 


^ sin 26 


( 60 ) 


The parameter X is evaluated by minimizing the modal strain 
energy. 

In the rigid body mode and the shear beam mode, each 
cross-section (z = constant) moves as a rigid body in the 
athwartship direction. In the buckling mode a cross-section 
is deformed, but the tangential strain = 0 . 

A. ADDED MASS EFFECTS 

Conceptually the effects of fluid pressure p^^ resulting 
from the radial motion of the structure should be assigned 
to the generalized force vector. It will be shown that the 
corresponding components are linear functions of q^, ^'2 and 
q^ so that it is appropriate to Include these effects by 
augmentation of the mass matrix. 
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Consider first a shell segment of infinitesmal length 
(dz) for which the radial displacement is 


u = cos e + U 2 cos 29 , 

(61) 

where = -q^ - q 2 sin ^ , 

(62) 

, . 2 irz 

and U 2 = -q^ sin — 

(63) 

The pressure p exerted on the hull as a result 
m 

of 

radial motion is determined using potential theory. 

For an 


incompressible fluid 


< 

B 

II 

0 

• 

(611) 

Let Pm “ Pi P 2 • 

(65) 

If Pi “ Pq f ® ’ 

(66) 

2 

and P 2 ~ Pq ^ ’ 

r 

(67) 

then p^ satisfies equation (6^]). 


'From Euler's equation of motion for an invlscid 

fluid, 

•77— „ = - p u, cos 6 , 

9r r=a 1 

(68) 

p^ = p a u^ cos 9 , 

r=a 

(69) 
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where p Is the fluid density. Similarly, 


r=a 


= ^ cos 20 


(70) 


Therefore, the pressure is given by 


m 


r=a 


= [P;]L P2^r=a “ Pa(u^ cos 9 + ^ 

(71) 


Let ' be the work done by the pressure p^ during a 
virtual displacement 6u. 

Then 

2it 

6W ' = -dz / p 6u a d0 . (72) 

P 0 

Substituting 

6 u = 6 u^ cos 0 + 6 u 2 cos 20 , 


and performing the indicated integration gives 

p •• n ' 

6Wp' = - irpa (u^ 6u^ ^ 2^ ^^2 ^^2^ 

If the substitutions 

TT 2j 

6 u^ = - 6 q^ - 6 q 2 sin — , 

and 


jt f 4 2irz 

6 u 2 = - 6 q^ sin -rj— , 


(72) 


(73) 

(7^) 
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are made and the Integration over the length carried out 
we have 


«SWp = - [6q^,6q2,6q2] 


m 

iE 

IT 

2 m 

m 

IT 

2 

0 

0 


0 

0 




(75) 


2 

where m = pira I (the displaced mass). The added mass matrix 
is 


[Ma] 


m 

2 m 

IT 


m 

IT 

2 

0 

0 



(76) 


B. STRUCTURAL MASS MATRIX 

For a shell having mass y per unit circumferential length 
the kinetic energy in an axial length dz is taken to be 

n ^ • p • p 

T ' = 5 - dz / y(u^ + v"^) a d0 . (77) 

s ^ 0 

But, from equations (58) and (59), 

u = u^ cos 0 + U 2 cos 20 , ( 78 ) 

and V = -u^ sin 0 “ ’ (79) 


5^1 






























































where, from equations (62) and ( 63 )> 


• • • . 7TZ 

Ui = -q^ - q^ sin -5- 

U2 = -q3 sin^ ~ . 


(80) 

( 81 ) 


Making the appropriate substitutions and evaluating the 
resulting integrals over 6 and the length gives 


2T 


s 


m 



+ 




( 82 ) 


Where m = 2 iTyail 

= shell mass (including stiffening rings). 

Adding the bulkhead mass m^ to the rigid body mode and 
including the added-mass-matrix the total kinetic energy 
becomes 


2T 


= (m + m^ + m^) q^^ + ^ (m + m^) q^q^ 
+ I (m + m^) q^^ + ^ (12m + 15m^) 


( 83 ) 


For neutral buoyancy: m = m^ + m^^ . Let B = , then 

equation ( 83 ) can be written as 


2T = {q}”^ [M] {q} , 


(8i<) 


55 



























where 


[M] = m 






f (1+6) 


0 

0 


0 ^( 12 + 156 ) 


(85) 


C. STRUCTURAL STIFFNESS MATRIX 

The coefficients of the structural stiffness matrix are 
determined from energy considerations. The strain energy V 
is expressed in the form 

2V = {q}”^ [K]{q} . (86) 


Since there is no strain energy associated with the 
rigid body mode and the other modes are elastically uncoupled, 
the stiffness matrix takes the form 


(87) 


1. Shear Beam Mode 

The modal displacements are given by equations (59). 
There is no ring bending and shell bending contributions 
are negligible, so that only m.embrane energy need be 
considered, i.e., Cg = 0, 0 and 


[K] = 


0 

0 


k 


0 

22 
0 k 


0 


0 


33 
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rfij M 


(88) 


V = — + i — 

'z0 " 9z r 96 


IT . A 

j sin 6 cos -j- 


The strain energy associated with the shear beam mode is 

given by 


= i i « ‘'s * ■ 


(89) 


0 0 


or 


2V 


sb 


G hg a ir- 

2l 


(90) 


2. Buckling Mode 

The strain energy of the buckling mode is determined 
in two parts: the shell bending energy and the membrane 
energy The modal displacements are given by equations (60) 

The shell bending energy is given by 


2V^ = 


H 2ir 

/ / 

0 0 


(D„K„^ + 2vD„k Kq+DaKa^ + ^D aKa„^) h_a d0dz 
zz zz6 66 z66zs 

(91) 


vfhere and D. are the flexural rigidities given as 
z u 


E h 


D = 


(12(1 - v'^) 


(92) 


and 


^6 ^ ^^z • 


(93) 
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The value a is a constant of proportionality determined by 
considering the orthotropic properties of the ring reinforced 
structure. The quantity is the 'twisting' rigidity 
given by 


D 


Qz 


1 - V 


D, 


. (9^1) 


where effects of ring twisting are neglected. 


The curvatures ic^, Kq and are determined from equations 


(60) as: 


K = ^ 

" 3z2 


*=6 " 


[ ^ + u ] 


a ae‘ 


and 


(95) 


2 

1 8 u 


ez a 303z 


Substituting and performing the indicated integrations yields 


ttD 


2V^ = 


b 2 


12-n^ ^ ^ a i + H(i - v)Tr2 f: + 


I I ^3 


(96) 


The membrane energy is 


1 217 


2V, 


m 


= / / (E'e^^ + Gy^e^) h^a dSdz , 


0 0 


(97) 
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v;he re 


e 

z 


3v; 

9z 


2ttX 
i ^3 


cos 


2Trz 

a 


cos 


26 


(98) 


Y 


ze 


3v ^ ^ 

3z r 36 


^3 ^ 2il 


2X. 


sin 


2ttz 

£ 


sin 


26 


(99) 


E' = E/(l - v^) 


Performing the indicated substitutions and evaluating the 
integrals gives 


2V 


m 


irailhg 

2 


IlTr2x2 

il2 


E' 


+ 


JL . 2X 
a 


G] q32 


( 100 ) 


The parameter X is determined by minimizing the 


membrane energy with respect to X, i.e., 
result of this minimization is 



The 


X 


-ir 


£ 

a 


^ (ir‘ 


E' 

G 



( 101 ) 


The required coefficients of the stiffness matrix [K] 
are determined from equations (90), (96), (100) and (101). 

D. GENERALIZED FORCES 

The generalized force vector {f} is determined by employing 
the virtual v-jork principle. Let 6W be the work done by the 
pressure p during a virtual displacement 6q. Then 
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2n I 

6^ = - f f p (|^ 6q 
n n ‘’Ht 


. 9u x.„ 
1 9qo 


9u 


2 9qo "'"S 


6q_) a dzde (102) 


v/here 


9u 


3q^ ’ 9q2 


and 


9q. 


are determined from the modal 


displacement equations (58-60), i,e,, - 


9u 

9qn 


= cos 0, 


9U . TTZ _ 9u 

- 35^ ° X ® ’ - 35^ 

the integral over the surface 


= sin^ ^ cos 20 


Evaluating 


.6W = [6q^,6q2,6q23 





where 


f 


f 


f 


1 

2 

3 


2ir 


aS, f p 


0 


2ail 

TT 



2it 

/ 

0 

2t\ 

P 


cos 0 d0 , 


p cos 0 d0 , 


cos 20 d0 


(103) 


(10i|) 


(105) 

( 106 ) 


The generalized forces are calculated from the two- 
dimensional finite element pressure-time history. The 
resulting genex’alized forces are given in Fig. 1^1. The 
rigid body force f^ is in close agreement with the results 
calculated by Sette [1] and Murray [3] for a step pulse. 
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Pig. 1^1. Generalized Forces-v.s .-Ti.'iie 
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E. GENERALIZED DISPLACEMENTS 


The Newmark-6 method, equation (42), of time integration 
was employed to evaluate the time dependent response of the 
structure. Figure 15 shows the resulting generalized 
displacements (q^^^E/p^a) as a function of the nondlmensional 
time parameter ct/a. 

F. STRUCTURAL STRESSES 

The structural stresses are determined from the strain 
relations utilized to determine the coefficients of the 
structural stiffness matrix. 

For the shear beam mode equation (88) gives 


“ ■'ze “ °''ze = 0 1 qj sin e cos ^ 


(107) 


For the buckling mode the stress due to ring bending is 




. 2 TTZ 

q^ sin -j- cos 26 


( 108 ) 


v/here c^ is the distance from the centroid of the cross- 
section to the extreme fiber. The membrane tension stress 
(equation (98)) Is 


S 


3 


a 


z 


z 



2uz „» 

cos —^ cos 26 


(109) 
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and the membrane stress (equation (99) is 


“ ^z6 ~ ^ *^3 ^ 21 a ^ t ’ (110) 

where X is determined from equation (101). 

Table I shows the maximum nondlmenslonal stresses Sj^/p* 
as a function of ct/a. 



TABLE I. 

MAXIMUM 

MODAL STRESSES 


ct 

"l 

^2 



a 





0.33 

0. 049 

0.035 

-0.036 

0.042 

0*67 

0.720 

0.467 

-0.469 

0.553 

1.00 

3.242 

1.880 

-1.887 

2.223 

1.33 

8.561 

4.370 

-4.386 

5.167 

1.67 

16.644 

7.351 

-7.377 

8.691 

2.00 

26.700 

^ o 

X O • w “» c. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


The results of this investigation confirm that the finite 
element method provides a satisfactory numerical method of 
evaluating the interaction of a plane acoustic blast vjave 
and a rigid structure. Extension of the method to other tv/o 
and three-dimensional geometries using the techniques 
employed appears to be practical. 

The submarine structural model is a simple three degree- 
of-freedom model. It was chosen to provide an illustrative 
example of the implementation of the superposition theorem 
and not as a general model of a submarine. 

It is recognized that a more sophisticated model is 
necessary to determine the response of a subm.arine subjected 
to blast wave loading. To improve the model it would be 
advantageous to: 

a. Include additional modes, particularly a 'breathing' 
mode to incorporate the effects of hydrostatic loading; 

b. Distribute and elastically mount the mass of the sub¬ 
marine not Included in the shell; 

c. Remove the restrictions on the rigidity of the 
bulkheads. 

Incorporating the above improvements and nonlinear 
structural characteristic into the submarine model will 
provide for a more complete analysis of the dynamic response 
of a submarine to blast-wave loading. 
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APPENDIX A 


PRESSURE-vs.-TIME AT THE RIGID CYLINDER-FLUID 
INTERFACE FOR 6 = 30°, 60°, 90°, 120° and 150° 

The pressure-time histories at the rigid cylinder-fluid 
Interface at 9 = 30°, 60°, 90°, 120° and 150° are shown In 
Figs. 16-20. The solid line represents Sette’s [1] results 
displaced along the time axis to coincide with the incidence 
of the constant portion of the cubic ramp-step v;ave. 

Sette [1] provided no data for 0 = 120° and 6 = 150° and 
therefore no comparison is available. Data for 6=0° and 
9 = l80° is presented In Figs. 9 and 10. 

The integration, wave and mesh characteristics are: 

ch _ 1 
a F * 

c^r 

0.97 , 

In = i 

a 2 * 
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Fig. 16. Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (0 = 30®) 
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Fig. 17. Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (0 = 60°) 
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Fig. 18. 


Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (0 = 90°) 
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Pig. 19. Pressure-Ys.-Tir.e at the Rigid Cylinder- 
Fluid Interface (0 = 120°) 
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Fig. 20. Fressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (0 = 150°) 
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APPENDIX B 


TABULATED PRESSURE-TIME HISTORY 
AT THE RIGID CYLINDER-FLUID INTERFACE 


The following tabulated data is the result of a tv/o- 
dimensional finite element solution of the interaction of 
a ramp-step wave and a rigid cylindrical structure. 

The integration, wave and mesh characteristics are: 
h = a/6c. 


T 

r 


0 

1 

2 


97 a/c, 

a . 
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FLUID INTERFACE (continued) 























































p/p* AT INDICATED ANGLE AND TIME 
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TABLE II. PRESSURE AT THE RIGID CYLINDER- 
FLUID INTERFACE (concluded) 

































































APPENDIX C 


SUBMARINE STRUCTURAL PARAMETERS 

The specific structural parameters used in the calculation 
of the coefficients of equation (57) are: 

p = fluid density = 9.36 x 10~^ Ibf sec^/in^ 

E = 30 X 10^ psl, 

G = 11.6 X 10^ psi, 

V = 0.300, 

a = hull radius = l 80 in., 

H = bulkhead spacing = H80 in., 

h = shell thickness = 1.60 in., 
s 

$ = m /m = 0.200 

s 

The internal ring reinforcing frames are characterized 
by: 

Frame spacing = 36.1 in., . 

Web thickness = 1.10 in.. 

Frame depth = 9.20 in.. 

Flange width = 9.20 in.. 

Flange thickness = 1.80 in., 

a = orthotropic constant of proportionality = 127. 
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The mass matrix is 


[M] = 


9.^05 3.592 0. 

3.592 2.821 0 . 

0. 0. 1.102 


X 10^ Ibf sec^/in 


The stiffness matrix is 


[K] = 


0 . 

0 . 

0 . 


0 . 

107.9 

0 . 


0 . 

0 . 

2.236 


X 10° Ibf/in, 
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